% ConditionedMediaConcentrationProfile_3D.m
% Eden Ford
% 02 August 2025

% Code to model diffusion-reaction concentration profile into a cylindrical
% hydrogel

% Please excuse the lack of efficiency in this code. I am not a computer
% scientist.

close all
clear all
clc

SumLimitM = 100; % m describes columns
SumLimitN = 100; % n describes rows

delta_t = 1*60*60; % seconds; 1 hour time intervals for solution concentration update
C_i_max = 100; % mmol/cm3 (arbitrary concentration to represent percent CM); max outer concentration
writer = VideoWriter('CM_Video');
writer.FrameRate = 2;
open(writer)
%% Constants and variables
% All in cm, s

D_water = 9.38e-07; % cm^2/s; diffusivity in water for PEGF
D = 7.99e-07; % cm^2/s; effective diffusivity for PEGF
t_half_hr = 1.8; % hr; protein half life (PDGF)

R = 4.6/10; % cm; hydrogel radius
H = 2.0/10; % cm; hydrogel height

t_half = 2.4*60*60; % s; protein half life
k = log(2)/t_half; % 1/s; protein degradation rate constant
V = 0.25; % cm^3 (mL); solution volume

r_list = [0 : 0.005 : R]; % cm; range of r locations
z_list = [0 : 0.005 : H]; % cm; range of z locations

[r,z] = meshgrid(r_list,z_list);

%% Constant Vectors

EmptyMatrix = zeros(SumLimitN,SumLimitM);

% zeros of bessel function (sum over m)
j0_m = besselzero(0,SumLimitM,1)';

% mu (sum over m)
mu_m = EmptyMatrix;
for i = 1:SumLimitM
    mu_m(i,:) = j0_m(:)/R; % mu (sum over m)
end

% lambda (sum over n)
lambda_value = zeros(SumLimitN,1);
for i = 1:SumLimitN
    lambda_value(i) = (2*pi*i-pi)/(2*H); % lambda (sum over n)
end
lambda_n = EmptyMatrix;
for i = 1:SumLimitM
    lambda_n(:,i) = lambda_value(:);
end

%% Concentration: t=0
C_t0 = r*0+z*0;
t = 0;
C_inf = 100;
% Heatmap
figure
    hold on
    Plot1 = surf(r,z,C_t0);
    set(Plot1,'edgecolor','none')
    cmap = flip(gray);
    %colormap(cmap)
    %title('t = 0')
    xlabel('hydrogel radius [cm]','FontSize',12)
    ylabel('hydrogel height [cm]','FontSize',12)
    Bar = colorbar;
    ylabel(Bar,'conditioned media concentration [%]','Rotation',-270,'FontSize',12);
    hold off
    xlim([0 0.46])
    ylim([0 0.46])
    caxis([0 C_i_max])
    txt0 = ['Initial Conditions'];
    text(0.025,0.28,txt0,'FontSize',14,'FontWeight','bold')
    txt1 = ['Solution Concentration:'];
    text(0.025,0.25,txt1,'FontSize',14)
    txt2 = [num2str(C_inf) '% Conditioned Media'];
    text(0.025,0.22,txt2,'FontSize',14)
    frame = getframe(gcf);
    writeVideo(writer,frame)

%% Concentration Over Time
C_inf = 100; % mmol/cm3 (arbitrary concentration to represent percent CM); initial conditioned media protein concentration in solution
disp(['Solution Concentration: ' num2str(C_inf) ' % CM'])

% Initial constants
A = -(4*C_inf*D*(mu_m.^2+lambda_n.^2).*sin(lambda_n*H))./(H*R*mu_m.*lambda_n.*besselj(1,mu_m*R).*(D*mu_m.^2+D*lambda_n.^2+k)); % constant, A; units of concentration
B = (4*k*C_inf*sin(lambda_n*H))./(H*R*mu_m.*lambda_n.*besselj(1,mu_m*R).*(D*mu_m.^2+D*lambda_n.^2+k)); % constant, B; units of concentration

% Flux constants (doesn't change)
nT_A = (2*pi*R*D*besselj(1,mu_m*R).*sin(lambda_n*H).*(mu_m.^2+lambda_n.^2))./(mu_m.*lambda_n); % volume/s

PlotTime = [1 : 1 : 72]; % in hours
% Cycle through PlotTime
for b = 1:length(PlotTime)
    time = PlotTime(b); % plot timepoint in hours
    disp(['Timepoint: ' num2str(time) ' hours'])
    
    % Calculate concentrations at each time step, update C_inf and plot (each hour)
    t = delta_t;

    % (1) calculate time component
    TimeFunction = exp(-t.*(D*mu_m.^2+D*lambda_n.^2+k)); % unitless
    GammaFunction = A.*TimeFunction-B; % units of concentration
    
    % (2) calculate concentration at each location
    C_final = 0*z; % empty matrix
    for i = 1:length(z_list) % values along z
        for j = 1:length(r_list) % values along r
            z_i = z_list(i);
            r_j = r_list(j);
            C_rz = GammaFunction.*besselj(0,mu_m.*r_j).*cos(lambda_n.*z_i);
            C_rz_Sum = sum(sum(C_rz));
            C_final(i,j) = C_inf+C_rz_Sum;
        end
    end

    % (3) calculate "moles" CM through hydrogel surface over time delta_t
    TimeFunction2 = 1-exp(-t*(D*mu_m.^2+D*lambda_n.^2+k)); % unitless; time function
    InnerFunction = (A.*TimeFunction2./(D*mu_m.^2+D*lambda_n.^2+k))-B*t; % concentration*s
    n_T = nT_A.*InnerFunction; % mol; mol transport from solution over delta t
    n_T_Sum = sum(sum(n_T)); % mol; mol transport from solution over delta t
    
    % (4) update C_inf
    C_inf0 = C_inf;
    C_inf = (V*C_inf + n_T_Sum)/V; % updated C_inf, no solution degradation taken into account
    %C_inf = C_inf_k0*exp(-k*delta_t) % updated C_inf, with degradation taken into account
    disp(['Solution Concentration: ' num2str(C_inf) ' % CM'])
    C_inf_text = C_inf;
    
    % (5) update solution concentration at 24 and 48 hours
    if time == 24
        C_inf = 50; % update with fresh solution after 1 day
        disp('')
        disp('Replace Solution with Fresh Media')
        disp(['New Solution Concentration: ' num2str(C_inf) ' % CM'])
    elseif time == 48
        C_inf = 50; % update with fresh solution after 2 days
        disp('')
        disp('Replace Solution with Fresh Media')
        disp(['New Solution Concentration: ' num2str(C_inf) ' % CM'])
    end
    
    % (6) update constants
    Gamma0 = GammaFunction; % old C_inf
    B = (4*k*C_inf*sin(lambda_n*H))./(H*R*mu_m.*lambda_n.*besselj(1,mu_m*R).*(D*mu_m.^2+D*lambda_n.^2+k)); % new constant, B; units of concentration; new C_inf
    A = (((4*sin(lambda_n*H))./(H*R*mu_m.*lambda_n.*besselj(1,mu_m*R))).*(C_inf0-C_inf)) + B + Gamma0; % new constant, A; units of concentration; new C_inf
    
    % (7) display results and produce plots
    % Heatmap
    figure
        hold on
        Plot1 = surf(r,z,C_final);
        set(Plot1,'edgecolor','none')
        cmap = flip(gray);
        %colormap(cmap)
        %title(['t = ' num2str(time) ' hours'])
        xlabel('hydrogel radius [cm]','FontSize',12)
        ylabel('hydrogel height [cm]','FontSize',12)
        Bar = colorbar;
        ylabel(Bar,'conditioned media concentration [%]','Rotation',-270,'FontSize',12);
        hold off
        xlim([0 0.46])
        ylim([0 0.46])
        caxis([0 C_i_max])
        txt0 = [num2str(time) '-Hour Timepoint'];
        text(0.025,0.28,txt0,'FontSize',14,'FontWeight','bold')
        txt1 = ['Solution Concentration:'];
        text(0.025,0.25,txt1,'FontSize',14)
        txt2 = [num2str(C_inf_text) '% Conditioned Media'];
        text(0.025,0.22,txt2,'FontSize',14)
        frame = getframe(gcf);
        writeVideo(writer,frame)
    GelMinimum = min(min(C_final))
    GelMaximum = max(max(C_final))
    GelRange = GelMaximum-GelMinimum
end
close(writer)